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Hirsch’s dynamic Hubbard model describes the effect of orbital expansion with 
occupancy by coupling the doublon operator to an auxiliary boson. We use the 
nonequilibrium dynamical mean field method to study the properties of doublon 
and hole carriers in this model in the strongly correlated regime. In particular, 
we discuss how photodoping leads to doublon and hole populations with different 
effective temperatures, and we analyze the relaxation behavior as a function of the 
boson coupling and boson energy. In the polaronic regime, the nontrivial energy 
exchange between doublons, holes and bosons can result in a negative temperature 
distribution for the holes. 


I. INTRODUCTION 

Strongly correlated electron systems are often described in terms of the Hnbbard model 
or its mnltiorbital extensions. This model assnmes that the electrons interact only locally 
on each lattice site, and that this interaction is of a particle-hole symmetric form. In the 
case of the single-band Hubbard model, the interaction on site i, parametrized by U , can 
be written as Uijii^ — \){ni^ “ i)’ energy costs for adding and removing an 

electron in the half-hlled state are the same. As was pointed out by Hirsch and co-workes in 
a series of papers, this model ignores the fact that doubly occupied atomic orbitals expand 
as a result of Coulomb interactions and electron correlations. To capture this physics, 
Hirsch introduced the dynamic Hubbard model,- which contains an additional coupling of 
the doublon operator to an auxiliary boson. The expansion of the orbital and associated 
reduction of the Coulomb interaction is then described by the dynamics of the boson and the 
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resulting screening effect. Several variants of this model have subsequently been formulated 
and explored.-^- These equilibrium studies have revealed interesting properties resulting 
from the hlling-dependent particle-hole asymmetric correlation effects, which may even be 
relevant for the understanding of unconventional superconductivity.-'^ 

More recently, a dynamic Hubbard model has been proposed to explain the behavior of 
organic charge-transfer salts in which molecular orbitals are periodically modulated via the 
selective excitation of an optical phonon mode with intense laser helds.-i^ In this context, the 
particle-hole asymmetric electron-boson coupling provides an experimental “knob” to tune 
the interaction effects in the solid, and hence raises the interesting prospect of dynamical 
control of material properties. 

Theoretically, the nonequilibrium properties of dynamic Hubbard models are largely un¬ 
explored, but as in equilibrium, one can expect interesting effects resulting from the particle- 
hole asymmetric correlations. Examples are the nontrivial interplay between heat and par¬ 
ticle currents (Seebeck effect), or the unusual nature of the photo-doped insulator with its 
coexisting weakly/strongly correlated doublons/holes. Motivated by these considerations, 
we investigate here some nonequilibrium phenomena in the dynamic Hubbard model by 
means of the nonequilibrium dynamical mean held (DMFT) formalism.- In particular, we 
will study how the dynamic Coulomb correlations are rehected in the relaxation and recom¬ 
bination of doublons and holes produced by an impulsive stimulation of the Mott insulating 
state. We will show that the particle-hole asymmetry leads to different transient tempera¬ 
tures for the doublons and holes during the relaxation, and that the respective cooling rates 
depend in a nontrivial way on the coupling strength of the auxiliary boson. 

The rest of this paper is organized as follows. Section [TTl describes the model considered 
in this study and the approximate impurity solver used in the DMFT calculations. Sec¬ 
tion m presents results for the nonequilibrium spectral function and the time-evolution of 
the effective doublon and hole temperatures, while Sec. IIVI provides a short summary. 


II. MODEL AND METHOD 


We consider the dynamic Hubbard model introduced in the original paper by Hirsch 
^ Vijcl^Cja + ^ n^) + + V2gXini^ni^ , 


( 1 ) 
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where Vij denotes the hopping amplitude, U the onsite repulsion, /i the chemical potential 
and rii^ measures the occupation of site i with electrons of spin a. The coupling of the double 
occupancy to a harmonic oscillator with frequency cuq allows to describe the expansion of 
the atomic orbital (and associated reduction of the Coulomb interaction) which occurs when 
the orbital is occupied by two electrons. Xi and Pi are the boson position and momentum 
operators, respectively, and g is the electron-boson coupling. (Here, the expansion of the 
orbital would correspond to a shift of X in the negative direction.) In DMFT,— this lattice 
model is mapped to a self-consistent solution of a quantum impurity model with Hamiltonian 
H = Hioc + hfhyb + hfbath and a local term 


Hioc = Ud- p(nt + n^) + (a;o/2)(P" + X^) +V2gXd, 

'-V-" '-V-' 

= -^^electron =-^boson 


( 2 ) 


where d = n^ni measures the double occupation. The hybridization term i/hyb = 
+ h.c.) and the noninteracting electron bath i^bath = are de¬ 

termined by a self-consistent procedure in such a way that the bath mimics the lattice 
environment. In the action representation of the impurity model, the coupling to the envi¬ 
ronment appears as a hybridization term S'hyb = dtdt' cl(t)Ai^(t,t')c^(t'), and the 

hybridization function Acr(t,t') is determined by the DMFT self-consistency. 

To solve this impurity model we perform a unitary transformation^ which transforms 
operators O —?■ O = and thus shifts X ^ X — Xo. The transformed local 

Hamiltonian reads 


Hloc — .^electron + hfboson + ~ 2^0 ~ V^gdXf) + X{y/2gd — UJqXq), (3) 


where we denote transformed operators by a tilde. The electrons and bosons can be de¬ 
coupled by choosing Xq = so that the unitary transformation becomes O = W'^OW 

with 


W = 


ipV^d 
e “0 , 


( 4 ) 


and the terms involving Xq in Eq. (|3]) evaluate to ^Xq — ^/2gXod = —^d. The transformed 
interaction and chemical potential therefore read 

U = U-—,jl = iJ. (5) 

Uq 


dt) 


We also have to evaluate the transformed electron creation and annihilation operators, 
= = e The creation operator gives a nonzero result only if it 
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acts on a state with n„ = 0. In this case, the operator W on the right becomes 1, and the 

“0 


—iP ria 


operator on the left becomes e “o Hence, we find cj. = e “o and similarly 
= e Writing P = -^{b' — b) in terms of the boson creation and annihilation 

operators and separating the contributions from the occupied and empty a state, we finally 
get 


Cry = e “0 


+ a, 


( 6 ) 

(7) 


where we have introduced the Hubbard operators Ca'^ = UaCa^ and = {l—n^)ca'^. Using 
Eqs. (jSD and ([7]), the hybridization term of the impurity action can be written as ^hyb = 

+ (The 

remaining terms involving products C^C and C^C are not written explicitly because they 
do not contribute to the lowest order expansion considered below.) 

As introduced for the Holstein-Hubbard model in Ref. 1^, an approximation that captures 
the essential physics in the insulating phase and in the metal-insulator crossover regime at 
not too low temperature is to partially resum the diagrams of a hybridization expansion 
(non-crossing approximation, NCA),— and to decouple the electron and boson terms by 
factorizing products t c))(t)i?(t, t')co-(t') in the action and in the 

Green function, with the bosonic function (T is the contour ordering operator) 


Hbo 


- ,(t> - «<)^„) - cosh(/?c.„/2)) 


= —?exp 


( 8 ) 


Here, denotes the earlier/later time on the Keldysh contour. The same electron-boson 
decoupling has recently been shown to provide the optimal representation of an Anderson 
impurity model with dynamically screened interactions in terms of a model with static 
interactions.— In the case of the dynamic Hubbard model, we follow the same procedure 


and write Rhyb —t ClA(t,t')B(t,t')Ca- + ClA(t,t')Cfj 


Similarly, the contour- 


ordered Green function Ga(t,t') = —i{Tcry(t)cl^{t')) is factorized within the NGA like 

G.(t,t') = -*(rG.(t)Gi(t'))5(t',t) -*(ra(t)Gt(t')). (9) 


The correlation functions —i{TCry(t)Cl(t')) and —i{TC„{t)Cl{t')) are then evaluated using 
the standard NGA diagrammatic expressions,— as illustrated in Fig. [1] 
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FIG. 1: Illustration of the NCA self-energy diagrams for the dynamic Hubbard model. An 
open circle corresponds to an operator Ca = n^ct, a full circle to an open square to 

Ca = {l — na)ct and a full square to = (1 — na)C(j- A^- is the hybridization function, Q a 
pseudo-particle Green’s function, and B the bosonic propagator defined in Eq. ([8]). Red squares 
represent the “box weight” (I14p in a simulation with externally driven bosons. 


Because the dynamic Hubbard model is not particle-hole symmetric, an external driving 
of the bosons can be expected to yield nontrivial effects. Such a driving might have to be 
considered in the modeling of experiments similar to those recently reported by Kaiser et 
al.,— or Singla et al.,— where the lattice and hence the Wannier orbitals are shaken. This 
effect can be described within the present formalism, starting from the model 

.Hloc = ^^electron + ^^boson + {V^gd + f{t))X, ( 10 ) 


with f{t) some time-dependent function. Replacing V^gd —?■ '/2gd+ f{t) in Eq. ([3]), the 
Lang-Firsov shift becomes Xq = -^{\/2gd + and the shifts of the interaction and 

chemical potential are obtained from ^X^ — {\/2gd+ f)Xo = —^Xq = —i^{{g'^ + \/2gf)d+ 
^P)- Hence 

U{t) = U-\g^ + V2gf{t)), ( 11 ) 


and there is an additional time-dependent energy shift 


According to the discussion in Section HB of Ref. 
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, the time-dependent unitary trans¬ 


formation W{t) implies an additional term fH^(t)fkF(t) in the Hamiltonian. Including this 
term, the bosonic part of the transformed Hamiltonian becomes 


^boson = + ( 12 ) 

/ Uq 

From Eq. ([6]) and the solution of the Heisenberg equations of motion for the 6-operators, 
b^(t) = 6^6*“’“* -b one hnds (assuming /(O) = 0) 


( 13 ) 
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In the hybridization expansion, all vertices and C^. are thus multiplied by the “box 
weight” 

with s = ±1 for creation and annihilation operators (see Fig. d]). If the driving is of the 
form f{t) = y4sin(ci;f), then jldif{i) cos(a;o(t — i)) = (cos(ci;ot) — cos(a;t)) for oj ^ wq, 
and ^ sin(a;ot) for oj = ojq. 


III. RESULTS 


A. “Photo-doping” of the Mott insulator 


Since the approximate impurity solver described in the previous section is expected to 
give qualitatively correct results in the insulating phases of the model, we will consider a 
“photo-doping” type perturbation in a paramagnetic Mott insulator. For simplicity, we do 
not simulate the excitation of carriers across the Mott-gap by an electric held pulse, but 
instead apply a short U-pulse to the system: 


U{t) = 


U + AU sin(a;;7t) 0 < t < Att/oju, 


( 16 ) 


U t > Att/cuu- 

For uju = U, the effect of this modulation is to promptly create a certain number of doublon- 
hole pairs, so that the state of the system after the f/-pulse is similar to a photo-doped state 
after impulsive stimulation with a laser pulse with frequency comparable to U. We assume 
a semi-circular density of states of bandwidth 4u, so that the DMFT self-consistency loop 
simplihes to = v‘^Ga(t,t'),— and we use v as the unit of energy. 

The photo-doping induced changes in the spectral function of the (particle-hole symmet¬ 
ric) Holstein-Hubbard model have been investigated in Ref. 
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In contrast to the Hubbard 
model, where the Mott insulating spectral function is little affected by photo-doping,— in 
particular in the gap region, the rapid cooling of the photo-doped carriers and the pos¬ 
sibility of polaron formation in the Holstein-Hubbard model with strong electron-phonon 
coupling can lead to a series of doping-induced in-gap states. The appearance of these pola- 
ronic states also has an important effect on the life-time of the photo-carriers, which becomes 
much shorter than in the Hubbard model, and also doping dependent {d/dt{d) oc {{d) — dthY, 
with dth the doublon population of the initial state). 
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FIG. 2: Equilibrium spectral functions and compressibilities in the dynamic Hubbard model. The 
left panel shows spectral functions for U = 10, /x = 0, wq = 1 and indicated values of g, while the 
middle panel shows spectral functions for U = W, g = 0, fixed g'^ fuj^ = 1 and indicated values 
of loq. The right panel plots the compressibility as a function of g foic g = 1, loq = 1, U = 5 and 
U = 10. The inverse temperature is /3 = 5 for all plots. 


B. Equilibrium results 

We now turn to the dynamic Hubbard model ([I]) and first briefly discuss the spectral 
properties of the Mott insulating equilibrium system. The left panel of Fig. |2] plots the 
spectral function A{u)) for U = 10, /i = 0, cuq = 1 and indicated values of g. The inverse 
temperature is /3 = 5. Because the chemical potential is inside a large gap, the hlling is 
approximately n = 1 independent of g. For g = 0 (Hubbard model), the system is particle- 
hole symmetric and the spectral function features two Hubbard bands of approximately semi¬ 
circular shape. As the electron-boson coupling is increased, the spectral function becomes 
asymmetric: while the lower Hubbard band essentially retains its semi-circular shape, the 
upper Hubbard band broadens and for large enough g splits into sidebands. In addition, 
there is a small shift of the Hubbard bands to higher u. If the parameter g’^/too > 0 is kept 
hxed while Wq is varied, the separation between the sidebands changes, while the broadening 
of the upper Hubbard bands stays roughly the same (middle panel). 

The boson-induced broadening of the upper Hubbard band in the small -(7 regime is a 
manifestation of the reduced doublon correlations, which is the effect described in the original 
work by Hirsch. The splitting into sidebands at larger g is due to the formation of polarons. 
We will not discuss here under which conditions this polaronic regime may be realized, but 





























simply explore the nonequilibrium properties of the dynamic Hubbard model in the full 
range of boson couplings. 

One can also observe the asymmetry in the electron and hole correlations by computing 
the density n(/i) or the compressibility dn/dfi. As shown in the right panel of Fig. O the 
compressibility is larger on the electron-doped side than on the hole-doped side, which is 
consistent with weaker doublon correlations. 

C. Spectral functions of the photo-doped system 

We next investigate the changes in the spectral function after the impulsive excitation 
dehned in Eq. flT5]) . Figure [3] shows results for t/ = 10, = 1, Wq = 1 and pulse parameters 

Af/ = U and uju = U. The top left panel corresponds to a half-hlled system. We plot both 
the total spectral function A’"®* = -f- A^, the distribution of occupied states A^, and 

the distribution of unoccupied states A>. These time-dependent spectra were obtained by 
forward time integration, using the formulas 

1 r°° 

A<{u,t) = -lm (16) 

TT Jt 

A^^\u,t) = -Im / (17) 

TT Jt 

We plot the results for t = 15, i.e., during the relaxation which follows the short 2-cycle 
pulse (which lasts until An/cju = 1.26). As can be seen by comparing the transient spec¬ 
tral function to the initial equilibrium result in grey, the photo-doping induces substantial 
changes, especially in the gap region. In analogy to the previous results for the Holstein- 
Hubbard model^ we hnd that the production of doublons results in a series of in-gap states, 
which can be viewed as additional boson sidebands of the upper Hubbard band. However, 
in contrast to the particle-hole symmetric Holstein-Hubbard case, the sidebands are limited 
to the region a; > 0, and there are no in-gap states appearing in the vicinity of the lower 
Hubbard band. The lower Hubbard band merely exhibits an asymmetry reminiscent of a 
smeared-out quasi-particle band near the upper band edge. 

The distribution of occupied states shows that the additional doublons created by the 
pulse are rapidly cooled down and are the origin of the polaronic in-gap states. At the same 
time, the distribution of unoccupied states shows an accumulation of holes near the upper 
edge of the lower Hubbard band. The time-evolution of the doublon and hole populations 
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FIG. 3: Left and bottom panels: Nonequilibrium spectral functions measured at time f = 15 in a 
half-filled, hole-doped, and electron-doped system after photo-doping. The parameters are U = 10, 
5 = 1, (jQ = 1 and /3 = 5. The frequency of the two-cycle excitation pulse Eq. (fl^ is oju = U 
and the amplitude is AU = U. The black curve shows the total spectral function, while the green 
and pink curves indicate the distributions of occupied and empty states. For comparison, we also 
plot the initial equilibrium spectral functions. The top right panel shows the time evolution of the 
distribution functions for n = 1. 


for n = 1 is illustrated in the top right panel. The curves with colors ranging from green 
to blue show the function in the lower Hubbard band from t = 3 to 21, and the set 
of curves with colors ranging from red to dark-violette show the function in the upper 
Hubbard band for the same times. One notices that the doublon population exhibits almost 
no time dependence, which shows that the cooling of the doublons and the appearance and 
population of the polaronic in-gap states is so fast that it essentially takes place during the 
application of the pulse. This observation is consistent with the Holstein-Hubbard results in 


Ref. 


16l for the same parameter regime. The redistribution of the hole population is a much 
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slower process, which is expected because the holes do not couple directly to the bosons 
in the dynamic Hubbard model. Nevertheless, the cooling of the holes which is evident in 
this set of data is very fast compared to the relaxation in the Hubbard model {g = 0), 
where there would be almost no change in the energy distribution of doublons and holes 
on the timescales of this plot.— Apparently, the scattering with “cold” doublons provides a 
relatively efficient dissipation mechanism for the holes. In the next section, we will analyze 
this relaxation process more quantitatively. 

We can repeat similar calculations also for hole and electron doped systems. In the 
bottom panels of Fig. |3] we show results for n = 0.85 and n = 1.15. In the hole-doped 
system, the spectral function features a narrow quasi-particle peak at the Fermi energy. 
After the pulse, this peak broadens due to the effects of heating and photo-doping. Near 
the upper Hubbard band, we again observe the appearance of a series of polaronic side¬ 
bands. They are less prominent than in the half-hlled case, because in the doped system, 
the doublon-hole recombination is faster, so that the number of photo-carriers in the upper 
Hubbard band is substantially reduced by t = 15. In the electron-doped system, we observe 
a series of doping-induced side-bands already in the equilibrium state. After the pulse, 
the weight of these polaronic states is further increased, while in the lower Hubbard band 
we observe the appearance of a small “quasi-particle” peak. Again, because of the faster 
doublon-hole recombination, the changes in the lower Hubbard band are smaller than in the 
half-hlled case. 


D. Effective doublon and hole temperatures 

Except deep in the polaronic regime, the doublon and hole distribution functions quickly 
reach an approximately thermal shape in the energy region corresponding to the upper and 
lower Hubbard band. It is therefore instructive to analyze the relaxation process in terms of 
the effective “temperatures” of the doublons and holes. To measure these temperatures, we 
compute the nonequilibrium distribution function f{oj,t) = In equilib¬ 

rium, f{uj) corresponds to the Fermi function for the corresponding temperature. As shown 
in the top left panel of Fig. |4]for the model with g = 0.5 and cuq = 0.25, the photo-doped 
system does not have a Fermi like nonequilibrium distribution. Rather, there are two Fermi 
edges, one in the energy range of the upper Hubbard band (cu ~ 4 in the hgure), and one in 
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FIG. 4: Top left panel: Nonequilibrium distribution function = A^{io,t)/A^^^{uj,t) for the 

model with g = 0.5, ujq = 0.25 and indicated values of t. To illustrate the energy range covered by 
the Hubbard bands, we also plot A^°^{uj,t = 60) (black line). For the t = 12 curve, we show the 
fits to Fermi distribution functions by dashed lines (fitting range [-7,-2] and [2,7]). Top right 
panel: Effective electron and hole temperatures extracted from the Fermi function fits for g = 0.5, 
oiQ = 0.25. The dashed line indicates the temperature of the initial state, and the arrow the end 
of the 2-cycle pulse. Bottom panels: Effective doublon and hole temperatures for different boson 
couplings and boson frequencies. 


the energy range of the lower Hubbard band (cu ~ —4). By htting the distribution functions 
to a Fermi function in these energy regions, one can extract an effective doublon and hole 
temperature, and in principle also the corresponding effective chemical potentials. For the 
earliest time (t = 12) these hts are indicated in the hgure by dashed black lines. From these 
hts, one extracts = 12) 0.23 and = 12) 0.92. 

The time evolution of these effective electron and hole temperatures is illustrated in the 
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top right panel. Since the Fermi hts depend on the htting range, there is some nncertainty 
associated with the calculation of the effective temperatures, and the error bars in the hgure 
give a conservative estimate of this uncertainty. It is obvious that the doublons, which couple 
directly to the bosons, are rapidly cooled down to an effective temperature close to the initial 
equilibrium temperature of 1//3 = 0.2 (horizontal dashed line), while the holes are still in 
the process of dissipating their energy at the longest accessible time. We note that in our 
closed system, the energy injected by the pulse will lead to a thermalization at a temperature 
above the initial equilibrium temperature, so that the effective temperature of the doublons 
might increase on longer timescales. However, since the energy changes associated with the 
redistribution of spectral weight within the Hubbard bands and the cooling of doublons 
and holes during the initial stage of the relaxation can apparently be easily absorbed by the 
bosons, such an increase of the doublon temperature would be controlled by the doublon-hole 
recombination time, which is very long in large-gap insulators.— 

The bottom panels show similar results for different boson couplings and boson fre¬ 
quencies. Obviously, the cooling rates and the ratio between effective doublon and hole 
temperatures depends on these parameters in a nontrivial way. 

E. Dependence of the cooling rate on the electron-electron and electron-boson 

interaction 

While the coupling to the bosons provides a dissipation mechanism for the doublons, 
the cooling of the holes is an indirect process which involves doublon-hole scattering. Here, 
one could imagine scattering processes with virtual break-up of the doublon, d -|- h —)■ t 
+ i ^ d + h, oi the scattering between a high energy hole and a low-energy doublon, 
which transfers energy from the hole to the doublon without such virtual excitations. If the 
former processes are relevant, one should observe a dependence of the cooling rate on U, 
because the local energies of the states d + h and t + i differ by U. Figure |5] shows the 17- 
dependence of the effective doublon and hole temperatures for different boson couplings and 
boson energies. There is a small dependence of the effective hole temperature on U, which 
suggests a more rapid cooling in small-gap insulators, and thus a certain role of the hrst 
type of scattering processes in the initial stage of the relaxation. The effective temperatures 
at later times are however essentially independent of U. It thus appears that the dominant 
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FIG. 5: Effective doublon and hole temperatures in systems with photo-doping concentration 0.1, 
indicated boson couplings and boson energies, and different values of U. The upper curves show 
the hole temperatures and the lower curves the doublon temperatures. 


scattering mechanism is of the second type, i.e., energy transfer from hot holes to cold 
doublons without virtual excitations. 

We also analyzed the dependence on the parameters g and Uq and found that the cooling 
rate of the doublons and holes is essentially determined by A = /ujq. This is illustrated 
in the left panel of Fig. |6l which shows the effective doublon and hole temperatures for 
U = 15, photo-doping concentration 0.1, and different parameter choices corresponding to 
A = 0.25. The rate with which the effective temperature decreases is similar in the three 
cases, even though the boson frequencies and coupling strengths are very different, and the 
slowest boson oscillation period (~ 100) is longer than the plotted time range. 

For hxed Wq = 1, the cooling rate depends in a non-monotonous way on the coupling 
strength g (see middle panel). If g is small, the energy dissipation to the bosons becomes 
slow. As a result, the doublon temperature decreases slowly and the effective hole and 
doublon temperatures approach each other. In the intermediate coupling regime, the energy 
dissipation is most effective and both the doublon and hole temperatures decay quickly to 
a value close to the initial equilibrium temperature. In the polaronic regime, the doublons 
cool down very rapidly (although their effective temperature is no longer very well dehned 
due to deviations of /(w, t) from a Fermi distribution function, see Fig. [7]), while the transfer 
of kinetic energy from the holes to the doublons becomes increasingly inefficient, so that the 
holes remain hot. The inefficiency of the cooling via doubIon-hole scattering is explained 
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FIG. 6: Effective temperatures in systems with t/ = 15 and photo-doping concentration 0.1. Left 
panel: constant effective coupling strength A = jooo = 0.25. Middle panel: constant boson 
frequency wq = 1- The upper curves show the hole temperatures and the lower curves the doublon 
temperatures. Right panel: inverse effective temperatures for doublons and holes as a function of 
g. The dashed (full) curves show the result at t = 20 (t = 40). 


by the narrow width of the polaron sidebands, and the relatively large energy separation 
between them. 

The 5 f-dependence of the inverse effective temperatures is shown in the right panel of 
Fig. ini The most efficient cooling occurs around g ~ 0.6 (corresponding to A ~ 0.4). For 
smaller g, the doublon and hole temperatures increase, while for larger couplings, only the 
hole temperature shows a marked increase. Interestingly, aX g ^ 1.9, the holes reach an 
inhnite temperature (flat) distribution, and for even stronger couplings, a negative tem¬ 
perature distribution (population inversion). From the comparison of the dashed and solid 
lines, which show the results for t = 20 and f = 40, it furthermore follows that the holes 
“cool down” in this negative temperature state (|Teff| decreases), i.e., the inverted population 
becomes steeper with increasing time. 

Figure [7] shows the spectral functions and distribution functions for g = 1.5 and g = 2, 
corresponding to the polaronic regime with positive {g = 1.5) and negative {g = 2) effective 
hole temperature. By comparing the spectral functions for the two coupling strengths, one 
notices that in the g = 2 case, the doublons populate in-gap states which hll the entire 
gap region. Due to these in-gap states, the doublon-hole recombination rate is substantially 
enhanced compared to g = 1.5. A possible explanation for the inverted hole population is 
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FIG. 7: Spectral functions and distribution functions in the polaronic regime (U = 15, loq = 1, 
g = 2 and 1.5). In the right hand panels, we plot the distribution functions in the energy range 
corresponding to a sizable hole and doublon population. For the later time (t = 40), we also 
show the Fermi function fits from which the effective doublon and hole temperatures have been 
extracted. For g = 2, the holes exhibit an inverted population. The high frequency noise is due to 
the cutoff of the Fourier integral. 

that the energy in a recombination processes is partially transferred to the hole subsystem, 
where it is stuck because holes do no longer efficiently thermalize with the doublons in 
the polaronic regime. The recombination process is still entropically favorable (although 
the “entropy” of the holons decreases with increasing negative temperatures) if a sufficient 
amount of energy is transferred to the bosons via the doublons. The rather long-lived 
transient state with inverted hole population thus results from doublon-hole recombinations, 
and the complicated energy exchange between holes, doublons and bosons in the polaronic 
regime with large cuq. On longer timescales, we expect that the system will thermalize in a 
state with identical, positive, temperatures for the doublons and holes. 
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IV. SUMMARY AND OUTLOOK 

We have discussed a DMFT based method for the simulation of nonequilibrium phenom¬ 
ena in the dynamic Hubbard model.- This class of models has originally been introduced 
to capture the fact that atomic orbitals expand if occupied by two-electrons, which results 
in particle-hole asymmetric correlations. Dynamic Hubbard models have also recently ap¬ 
peared in the discussion of externally driven solids, where the orbital shapes are modulated 
periodically.- We have considered here the model originally discussed by Hirsch, where the 
particle-hole asymmetric correlations are generated by a term gdX, which couples the dou- 
blon number operator to the position operator of an auxiliary boson. 

One can think of numerous nonequilibrium set-ups in which the particle-hole asymmetry 
leads to interesting effects, which are absent in a conventional Hubbard model or Holstein- 
Hubbard model description. While the external driving of the boson in a homogeneous 
Holstein-Hubbard model would result in a trivial time-evolution, the driving of the boson 
in the dynamic Hubbard model may produce nontrivial effects, ranging from particle-hole 
excitations to dynamically induced boson-magnon couplings. If a temperature gradient is 
applied across a slab, the unequal doublon and hole mobilities should result in a nonzero 
particle current. 

In the present work we have focused our attention on photo-doped paramagnetic Mott 
insulators, where the asymmetric boson coupling leads to the rapid cooling of doublons and 
a hot distribution of holes, since the holes can dissipate their energy to the bosons only 
indirectly, via doublon-hole scattering. We have shown that a short time after a photo¬ 
doping excitation, the doublon and hole distributions in the upper and lower Hubbard band 
exhibit an approximately thermal shape, so that one can extract an effective doublon and 
hole temperature. We have analyzed the time evolution of these effective temperatures as a 
function of the model parameters. As in the Holstein-Hubbard model, the cooling rate seems 
to be determined by A = g'^/too, and essentially independent of U. The fastest decrease of 
the effective hole temperature is observed for intermediate coupling strength, while in the 
polaronic regime, the energy transfer from the hot holes to the trapped doublons becomes 
slow. If the energy spacing uq between the polaronic states is large enough, the photo-doping 
of the Mott insulator may even lead to a long-lived transient state with an inverted hole 
population. 
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A long-lived negative temperature state could lead to potentially interesting effects. In 
the weakly correlated metallic regime, it has been shown that a system with a dynamically 
generated population inversion behaves as if the Coulomb interaction had been switched 
from repulsive to attractive.— While the present situation with strongly correlated photo- 
doped holes is quite different, one might nevertheless expect different phases to occur for 
T > 0 and T < 0 when |T| becomes small. Whether or not the inverted hole population 
leads to symmetry breaking to unconventional ordered states is an interesting open issue. 

In the future, apart from studying driven systems and inhomogeneous set-ups, it will 
be interesting to extend the nonequilibrium DMFT simulations to other types of electron- 


boson couplings, such as the gdX^ coupling of Ref. Q. The extension of the Lang-Firsov 
scheme to these variants of the dynamic Hubbard model is not straight-forward and may 
require additonal approximations. An alternative route is the development of combined 
strong/weak-coupling impurity solvers along the lines of Ref. l20 
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